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Abstract 

Dynamical systems that exhibit diverse behaviors can rarely be completely understood using a single approach. However, 
[by identifying coherent structures in their state spaces, i.e., regions of uniform and simpler behavior, we could hope to 
study each of the structures separately and then form the understanding of the system as a whole. The method we 

£S| present in this paper uses trajectory averages of scalar functions on the state space to: (a) identify invariant sets in the 
; i state space, (b) form coherent structures by aggregating invariant sets that are similar across multiple spatial scales. 
ClnFirst, we construct the ergodic quotient, the object obtained by mapping trajectories to the space of trajectory averages 
of a function basis on the state space. Second, we endow the ergodic quotient with a metric structure that successfully 
[captures how similar the invariant sets are in the state space. Finally, we parametrize the ergodic quotient using intrinsic 
diffusion modes on it. By segmenting the ergodic quotient based on the diffusion modes, we extract coherent features 
in the state space of the dynamical system. The algorithm is validated by analyzing the Arnold-Bcltrami-Childress 

'^j*. ' flow, which was the test-bed for alternative approaches: the Ulam's approximation of the transfer operator and the 

Q computation of Lagrangian Coherent Structures. Furthermore, we explain how the method extends the Poincare map 
analysis for periodic flows. As a demonstration, we apply the method to a periodically-driven three-dimensional Hill's 
vortex flow, discovering unknown coherent structures in its state space. In the end, we discuss differences between the 
"j^ ergodic quotient and alternatives, propose a generalization to analysis of (quasi-)periodic structures, and lay out future 
^ research directions. 

^~ ~ 1 Keywords: coherent structures; dynamical systems; ergodic partition; diffusion modes; trajectory averages 
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1. Introduction 

Historically, a lot of work in dynamical systems fo- 
cused on understanding a single type of behavior in a 
model: elliptical zones, hyperbolic trajectories, chaotic at- 
tractors, mixing regions, etc. As computers become more 
powerful, and experimental methods, such as Particle Im- 
age/Tracking Velocimetry, more precise, we are able to 
study large, complicated systems, whose state spaces com- 
prise many different coexisting coherent structures: re- 
gions of uniform and simpler behavior, which might be well 
understood individually. To apply our understanding to 
the entire system, however, we first need to identify where 
coherent structures lie in the state space. This task can 
be made difficult by the quantity of coherent structures, 
fractal boundaries between them, or the high dimension of 
the state space. 

Over the past decade, several different interpretations 
of coherent structures have appeared, each serving as a 
foundation for algorithms which extract them from the 
flow or from trajectory data. Coherent structures have 



been identified as: (i) distinguished trajectories that orga- 
nize behavior of the flow nearby, (ii) regions of state space 
between barriers of dynamical transport, (m) sets of initial 
conditions that are not dispersed by the flow, (iv) sets of 
trajectories on which all flow-invariant functions take con- 
stant values. All the mentioned definitions identify sets 
that arc dynamically invariant, in some sense, however, 
the choice of atomic objects, i.e., barriers of transport vs. 
initial conditions vs. trajectories, influences the type of 
approximation we obtain when using algorithms finite in 
time, space, and precision. 

Identification of coherent structures based on barriers 
to transport includes both classical geometric studies of 
invariant manifolds and newer approaches that attempt 
to generalize invariant manifolds to aperiodic, transient 
flows. The study of Lagrangian Coherent Structures (LCS) 
(lil . 17 1 focuses on structures that take the roles of hyper- 
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bolic invariant manifolds in time-varying flows. As a con- 
sequence, they are dominantly studied through finite-time 
Lyapunov exponent fields whose local extrema are tracked 
to obtain proxies to LCS [H, HIL H3 ■ 

While LCS-based methods identify distinguished hy- 
perbolic trajectories as organizing structures, a recent pa- 
per [28j . coauthored by the second author, proposed to 
identify both hyperbolic and elliptic behavior in a finite- 
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time, transient context. The approach, termed mcsochronic 
velocity analysis, studies the Jacobian of the flow averaged 
over short segments of trajectories. 

Both LCS and mesochronic velocity fields identify dis- 
tinguished trajectories that organize surrounding behav- 
ior. "Surrounding" is usually interpreted in terms of re- 
gions between neighboring distinguished trajectories. Such 
an interpretation is difficult to generalize to state spaces of 
higher dimension, where trajectories are not of codimension- 
one and cannot be used to segment the state space, with- 
out additional assumptions, e.g., symmetries of solutions. 

Approaches that aggregate smaller flow features into 
larger coherent structures are the alternative to analyzing 
barriers that separate trajectories. A prominent aggrega- 
tion method is the computation of invariant sets through 
cigenfunctions of the transfer (Perron- Frobenius) operator, 
a linear operator that captures evolutions of densities of 
points carried by dynamics [lol4l3j ). Numerically, the ap- 
proach most commonly proceeds by the Ulam's method: 
approximation of the system by a stochastic Markov chain, 
acting on a discretized state space. The eigenvectors of the 
Markov chain transition matrix approximate the eigen- 
functions of the transfer operator. The Ulam's method 
requires a computation of short trajectories started from 
a large number of initial conditions within each discretiza- 
tion cell. Despite its prominence, Ulam's method can be 
difficult to apply to, e.g., high-dimensional systems, and 
some experimental setups. A large number of state space 
dimensions may render computations infeasible, since the 
number of elements in Ulam's matrix grows exponentially 
with the dimension of the state space. On the other hand, 
if resetting and precise preparation of initial conditions is 
not possible, e.g., in certain fluids experiments, it can be 
more feasible to run fewer experiments for longer periods 
of time. 

The method that does use longtrajectories as a direct 
input has been pursued in 271 bill ■ It groups trajectories 
based on values of invariant functions along them. This 
approach interprets invariant functions as eigenfunctions 
of the Koopman operator, which captures how values of 
functions evolve along trajectories of the system. Koop- 
man and the aforementioned Pcrron-Frobenious operator 
form a dual pair in the appropriate function spaces. Know- 
ing eigenfunctions of the Koopman operator allows us to 
detect invariant sets in the state space, as level sets of 
Koopman eigenfunctions divide the state space into in- 
variant partitions. It is possible to compute them without 
approximating the operator itself: starting from any ob- 
servable, i.e., a scalar function on the state space, and 
averaging its values along trajectories projects the observ- 
able to the invariant eigenspace of the Koopman opera- 
tor. Choosing observables from a function basis, averaging 
them, and forming joint level sets of averaged observables 
results in finer and finer stationary partitions of the state 
space. Ultimately, the process converges to the ergodic 
partition of the space: the unique, finest partition into in- 
variant stationary sets. It was noted in (25j that even a 



small number of observables can be used to approximate 
the ergodic partition well. It remained unclear, however, 
how to select such a small set of observables that captures 
a lot of detail except by trial and error. 

The algorithm proposed in this paper follows the anal- 
ysis of dynamics using averaged observables. It docs not 
aim to "guess" an observable whose time average reveals 
the most information about the coherent structure; rather, 
it constructs a suitable set of invariant functions from av- 
erages of a function basis. To explain the method, we ex- 
pand on the representation of dynamics using the ergodic 
quotient [25| : each trajectory is represented by a sequence 
of time-averaged observables taken from a Fourier func- 
tion basis. Collection of all such sequences, the ergodic 
quotient, is equivalent to the ergodic partition, but more 
practical to work with, as analytic tools on sequence spaces 
are very well studied. To analyze the ergodic quotient, we 
endow it with a distance function, adapted from [26j , that 
identifies two trajectories as coherent if they spend simi- 
lar fractions of their evolutions within any particular set 
in the state space. As a result, the problem of dynamical 
analysis of flows is converted into a geometric analysis of 
the ergodic quotient, with a metric structure given by the 
empirical distance. Our approach is similar in spirit to 
analysis of integrable systems through geometry of their 
integrals of motion Q. 

To analyze geometry of the ergodic set we used the 
Diffusion Maps algorithm [fjj], which retrieves dominant 
modes of diffusion on objects with metric structure, e.g., 
branched manifolds. Wc introduced the use of Diffusion 
Maps as a tool in 0] ; here we explain its results in terms of 
the ergodic quotient. To validate the algorithm, we apply 
it to the Arnold-Bcltrami-Childrcss (ABC) flow, which was 
previously used as a test-bed for algorithms based on the 
Ulam's method [l2[ and on LCS jTij]. The second example, 
the periodically-forced Hill's vortex flow, demonstrates the 
application of the algorithm to time-dependent flows. 

When applied to analysis of time-periodic dynamical 
systems, our method extends the familiar analysis of the 
Poincarc map, including the averaged information about 
trajectories between two piercings of the Poincarc surface. 
Building on the application to time-dependent systems, 
we explain how the algorithm could be easily extended to 
the analysis of invariant features of periodic and quasi- 
periodic sets in the state space 27j. In this case, instead 
of ergodic averages along trajectories, we use harmonic 
averages along trajectories, which can also be interpreted 
through eigenspaces of the Koopman operator. While 
the ergodic averaging projected functions onto invariant 
eigenspace of the Koopman operator, the harmonic aver- 
aging projects functions to eigenspaces associated to other 
eigenvalues. 

There are several benefits to the analysis based on aver- 
aging observables along trajectories. Since we use the en- 
tire trajectory evolution, we do not need to access the en- 
tire state space directly, e.g., we can place initial conditions 
on a lower-dimensional surface and still obtain the full in- 
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formation about coherent structures intersected by that 
surface. Furthermore, compared to Ulam's method, we re- 
quire fewer initial conditions for analysis. These features 
makes the method applicable to experimental setups where 
resetting and initialization of experiments are restricted 
or costly. Moreover, when there are low-dimensional sur- 
faces of interest, such restricted initialization of trajecto- 
ries makes it possible to analyze even high-dimensional 
systems without the exponential growth in number of ini- 
tial conditions. Finally, approximation of eigcnfunctions of 
the Koopman operator does not require passing through 
a stochastic approximation of the deterministic dynamical 
system; instead projections of observables to eigenfunc- 
tions are evaluated directly by averaging. 

Although both of our examples in this paper origi- 
nated as models of fluid flows, the ergodic quotient analysis 
can be applied to a broad range of dynamical systems on 
compact state spaces. The incompressible fluid flows and 
hamiltonian systems, however, are two classes of systems 
that commonly have many invariant sets in their state 
spaces, arranged in intricate structure, therefore, those 
classes provide the most natural targets for our analysis. 
Nevertheless, dissipative systems can also be studied using 
the ergodic quotient, in which case the method identifies 
how basins of attraction are arranged in state spaces. 

The paper is structured as follows: the Section[2]presents 
the theory of the ergodic quotient, the empirical distance, 
and the diffusion coordinates, followed by the application 
to the autonomous ABC flow. Section [3] demonstrates 
how application to non-autonomous, periodically forced 
flows improves on the classical analysis of the Poincare 
map and supports the case by a numerical analysis of the 
periodically-forced 3D Hill's vortex flow. In Section [4] we 
discuss how our method compares with alternatives, as 
well as present our suggestions for application and future 
research. The Appendices analyze numerical properties of 
the averaging algorithm and the empirical distance, and 
present details about the Diffusion Maps that are neces- 
sary for an implementation as a computer code. 

2. The Ergodic Quotient 

Geometry of the ergodic quotient enables comparison 
of trajectories of dynamical systems across multiple spa- 
tial scales. To be able to study it, we move through three 
different spaces: the state space, the space of averaged ob- 
servables, where the ergodic quotient is constructed, and 
the space of diffusion coordinates, where the ergodic quo- 
tient is described by independent, intrinsic coordinates. 
Before delving deeper, we give an overview of the process, 
with italicized terms defined later in the text. 

To measure how similar two trajectories are in the state 
space (Figure [T]a), we use the empirical distance, which 
compares trajectories based on their residence times in 
spherical sets in the state space. On the state space, how- 
ever, the empirical distance is not practically computable, 
so we map trajectories to a space where it is. By selecting 
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Figure 1: A sketch of (a) the state space portrait, (b) the 
ergodic quotient in the space of averaged observables, and 
(c) the ergodic quotient in the space of diffusion coordi- 
nates. Crosses depict initial conditions used for numerical 
computation. Note that crosses on the ergodic quotient 
represent entire trajectories started from associated initial 
conditions on the state space. 

observables fk from a Fourier function basis and averag- 
ing them along trajectories, each trajectory gets mapped 
to a point on the ergodic quotient, which comprises se- 
quences of trajectory average^ fk for every trajectory in 
the state space (Figure[T]b). We endow the sequence space 
with a metric induced by a H~ s Sobolev space norm. The 
H~ s metric is topologically equivalent to the empirical dis- 
tance, yet it has a benefit that it is easily computable 
as a weighted euclidean metric. The ergodic quotient is 
typically low-dimensional, therefore its representation in a 
(infinitely-dimensional) sequence space is not economical. 
To address this issue, we use a nonlinear change of coor- 
dinates to represent the quotient in diffusion coordinates, 
preserving the intrinsic geometry of the quotient. Diffusion 
coordinates are an orthogonal, efficient, scale-ordered, low- 
dimensional representation of the ergodic quotient that 
has an effect of straightening the ergodic quotient (Figure 
[T]c). The euclidean distance in diffusion coordinates rep- 
resents the diffusion distance along the ergodic quotient. 
Diffusion distance is an intrinsic distance, locally consis- 
tent with empirical distance; for the sake of introduction, 
it can be thought of as a robust version of the geodesic 
distance along the quotient. 

Computationally, the construction starts with a finite 
set of trajectories Xi on the state space. They are mapped 
to the ergodic quotient by averaging a finite subset of the 
Fourier basis along them, therefore, the points on the nu- 
merical ergodic quotient are vectors of trajectory averages. 
Based on trajectory averages, we compute the matrix of 
pairwise H~ s distances between trajectories, which forms 
the input data to the Diffusion Maps algorithm. It com- 
putes the diffusion coordinates of each trajectory by solv- 
ing a matrix eigenproblem. Starting trajectories are then 
grouped into coherent structures, using a fc-means algo- 
rithm applied to the diffusion coordinates of trajectories. 
To visualize the coherent structures in the state portrait, 



Terms time average and trajectory average are used interchange- 
ably. In literature, ergodic average and Cesaro average are also used. 
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we either color trajectories based on their A:-cluster mem- 
bership, or we color them based on a single diffusion coor- 
dinate. 

2.1. Construction: From state space portrait to the ergodic 
quotient 

Take a compact euclidean manifold M. £ M. D as the 
state space of dynamics. The flow map is a composition 
(semi)-group of continuous functions (ft(x) : M. — > A4, 
with time t 6 T. For ODE-generated flows, the time-like 
parameter is a real number T = Rq , and the semigroup 
is assumed to be continuous in t. In map-generated flows, 
time progression over T = is interpreted as the iterated 
application of the map <p(x), i.e., ifit(x) = tp[ipt~i(x)], with 
(fo(x) = x. The space of observables J 7 is a set of Borel- 
measurable functions / : M. — > C; in this paper we will 
take T as the set of continuous functions on M . 

The geometric analysis of dynamics studies the state 
space portrait of the system, i.e., the collection of all tra- 
jectories in the state space, {(<Pt(x)) teVi ■ Va; £ M}. Chaotic 
behavior, however, implies that individual trajectories are 
non-robust to errors in initial conditions x. Despite the 
lack of robustness of trajectories, the averages of observ- 
ables along them are robust, even in chaotic regimes [1,|34[. 
As our approach is based on averaging, we expect the com- 
putations to remain robust even in chaotic regimes, which 
is the main motivator for this approach. Presented intro- 
ductory topics follow the exposition given in [22|, §4] . 

A finite-time average of an observable / is given by 



[/ o tp t ] (x)dt. 



Let E C M. be the set of points on which the limit as 

T — > oo exists. On E, we can define the infinite time 
average 

f(x) = lim - / [/ o ipt] {x)dt. 

If the flow map preserves a measure A, the set A^\E will 
be negligible in measure, A(.M\E) = 0, by Birkhoff's Er- 
godic Theorem. As hamiltonian systems and divergence- 
free flows preserve the Lebesgue measure, it is unlikely 
that an initial condition for which time averages do not 
converge will be selected when analyzing such systems. 
The same argument can be made for systems that pre- 
serve a Lebesgue-continuous measure, which contains an 
even larger class of systems. Consequently, we will equate 
E and Ai for purposes of this paper. 

The infinite-time average / i— > f(x) is a positive, lin- 
ear, continuous functional, and the Riesz representation 
theorem asserts that the averaging functional can be rep- 
resented by a Borel probability measure fj, x , the empirical 
measure, attached to the initial condition x: 



The empirical measures of sets fi x (E) can be interpreted 
as the fraction of evolution that the trajectory evolving 
from x spends inside the set E (residence or mean-sojourn 
time). For regular orbits, empirical measures are masses 
supported on fixed points or spread out along periodic or- 
bits, while on chaotic sets they may have a positive-area 
support. It is worth noting that empirical measures \x x are 
instantaneous invariants of the system, despite our use of 
an asymptotic process in their construction. They play a 
similar role in the measure-theoretic description of dynam- 
ics as the trajectories do in the geometric description. 

Instead of working with potentially singular empirical 
measures directly, we will work with their weak represen- 
tatives. Fix a countable basis of continuous functions fk 
on M., indexed by a multi-index fc, e.g., Fourier harmonics 
with D-dimensional wavevectors k £ Z D . Let the ergodic 
quotient map be it : M — > l°°(Ti D ), given by 

n(x) :=(..., f k (x),...) k , keZ D . 

Action of 7r at x can be interpreted as forming a weak 
representation of the empirical measure [i x . The set of all 
such sequences, generated by averaging from each initial 
condition, is the ergodic quotient £ := ir(M.). It will play 
the same role that the state-space portrait plays in geo- 
metric analysis, with the added benefit of robustness in all 
dynamical regimes. 

Strictly, every point on the ergodic quotient should be 
thought of as the set of trajectory averages of all contin- 
uous observables: by selecting a basis we merely pick a 
representation of the ergodic quotient in the particular se- 
quence space l°°(Z D ). It was shown in [27} that all such 
representations are equivalent; therefore, we will speak of 
any of them as the ergodic quotient when there is no am- 
biguity. 

Conceptually, we have moved from points in the state 
space x to associated averages of observables in two steps: 



1 11 / ; i \ 
x -> fi x — » (. . .,fk{x), . 



)k£Z D - 



(2) 



M 



fdfi x = f(x). 



(1) 



The representation steps (I) and (II) in © are not bijective 
in general. 

Step (I), at the very least, associates all points on the 
same trajectory with the same empirical measure; how- 
ever, the pre-images of empirical measures might contain 
more than a single trajectory. In Q we discussed the ap- 
proximation of the ergodic partition, which is the (unique) 
partition of the state space into ergodic sets. An invariant 
subset of the state space S C M. is called ergodic when the 
flow (p t restricted to it, ipt ■ S H> S, is an ergodic system. 
Since trajectory averages of ergodic systems do not depend 
on initial conditions, it follows that all points x £ S have 
the same empirical measure. In this sense, the pre-image 
of the ergodic quotient is the ergodic partition of the state 
space. Such an outcome is desirable, as ergodic sets are 
minimal robust atomic objects that contain (non-robust) 
trajectories. 
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Step (II) is bijective only if the entire function basis 
can be used to construct the ergodic quotient. On contin- 
uous state spaces which have infinite bases, we will need to 
truncate the function set to implement the algorithm on 
a computer: the magnitude of error introduced by basis 
truncation will depend on the distance structure on the 
ergodic quotient, discussed in the Section l2~2l 

Remark. Unless we are looking for the ergodic quo- 
tient of the entire system, we can place initial conditions 
only in a region of space that is of a particular interest. 
Moreover, only one initial condition from an ergodic set 
is required to represent the empirical measure for the en- 
tire ergodic set, as any trajectory in the ergodic set will 
correctly sample the associated empirical measure. There- 
fore, it is sufficient to sample the initial conditions from 
the surface that intersects coherent features we wish to ex- 
plore. After selecting the initial conditions and the set of 
obscrvablcs to average, we compute trajectories starting 
with each initial condition until time averages along them 
converge (see | Appendix A| , thus obtaining the approxi- 
mate mapping of the state space into the ergodic quotient. 

2.2. Metric structure: Empirical distance 

Trajectories are aggregated into coherent structures by 
a criterion of similarity embodied in a distance function 
between trajectories. The choice of the distance function 
determines the type of coherent structures we obtain. In 
this work, we chose to work with the empirical distance, 
which compares two trajectories based on residence times 
(2a |. Instead of working with the empirical distance di- 
rectly on the state space, we use a distance on the ergodic 
quotient that is equivalent to it, thus converting the anal- 
ysis of dynamics into analysis of geometry of the ergodic 
quotient. 

The empirical distance was originally used to quantify 
how closely a trajectory samples an a priori distribution on 
the state space, i.e., to compare an empirical measure with 
a fixed prior measure on the state space. The same dis- 
tance can be used to compare two empirical measures and 
estimate how similar the processes that generate measures 
are. Define the empirical distance V between measures [i x , 

by 



be used instead. The distance 



o J M 



{fe [B{p, r)] - [i v [B(p, r)]} 2 dpdr 



where B(p, r) C M. are spherical sets of center p and radius 
r normalized such that B(p, 1) D M. and Vol Ai = 1. Note 
that integration is over all spherical sets, without normal- 
ization by the sets' volumes, resulting in higher sensitivity 
to differences across larger spatial scales. 

Integration over uncountably many spherical sets is 
practically infcasiblc, however, the metric induced by the 
norm on a negative-order Sobolev space of measures can 



P_ s (/Xx,My)= F(x)-F(y) 



2,-s 



F(x) - F(y) 



E 



\fk(x) - h(y) 



(3) 



is easily computed as a weighted sum of Fourier coeffi- 
cients, F(x) = ( . . . , fk{x), • • • ) of measures involved. 

The order s of the Sobolev space is selected based on the di- 
mension D of the domain of measures Ai as s = (D + l)/2. 
Using 23_ s instead of T> does not change the resulting 
topology, as the empirical metric and the H~~ s with the 
order s = (D + l)/2 are equivalent, i.e., there exist posi- 
tive constants c and C such that for all [i x and fi y it holds 
that 

When fi x and [i v are empirical measures, the Fourier 
coefficients can be computed as trajectory averages of ob- 
servables, due to (JT|), which is the core observation on 
which we base our method. Therefore, we choose observ- 
ables fk as Fourier basis functions, 



hix) = (2tt) 



-D/2 



exp 



D 

i2TT^k d x d 

(2=1 



(4) 



with kd, and Xd being components of the wave- and state- 
vectors, respectively, assuming Ai = T D for simplicity 0. 
Interpreted in the ergodic quotient framework, sequences 
fk(x) and fk(y) are points on the ergodic quotient corre- 
sponding to empirical measures [i x and /i y . 

For purposes of our method, it is important to capture 
the topology and the geometry of the ergodic quotient, 
not the values of involved quantities. Therefore, we will 
use the empirical distance T> and the H~ s Sobolev space 
distance 2?_ s interchangeably: the former in explaining 
conceptual implications of the metric structure, the latter 
in computational considerations. 

In practice, we cannot compute averages of the entire 
function basis. The direct effect of the basis truncation is 
difficult to quantify; intuitively, the higher the cutoff K is 
set, restricting wavevectors to k E [-K, K] D , the finer are 
the state space features we are able to resolve through em- 
pirical distance. In |Appendix B| we show that truncation 
introduces at most 0(l/y/K) error in the distance for any 
state space dimension D. 

Remark. From a signal processing perspective, the 
H~ s distances correspond to a low-pass filtering of em- 
pirical distributions prior to their comparison. Since nu- 
merical errors are more likely to manifest in averaging of 



2 For arbitrary domains, observables are solutions of the Helmholtz 
problem on the domain. G. Mathew and I. Mezic will address the 
general formula for the H~ s distance in their upcoming work. 



5 



highcr-wavenumber modes, we expect that the use of such 
distances will improve the numerical stability of the anal- 
ysis. For a brief analysis of these effects, see | Appendix A| 

2.3. Geometry of the ergodic quotient 

Equipping the ergodic quotient with a distance func- 
tion provides a setting for studying the geometry of the 
space of invariant functions. The study of dynamics through 
ergodic quotient shows similarities with the study of dy- 
namics of intcgrable systems through geometry of integrals 
of motion. The main difference is that the ergodic quo- 
tient is a constructive technique: we do not need to know 
the expressions for integrals of motion to construct the er- 
godic quotient or approximate it numerically. Moreover, 
construction of the ergodic quotient is possible even when 
the system is not integrable. 




Figure 2: Geometry of the ergodic quotient under empiri- 
cal distance. Sketches of matching features in state space 
portraits (top) and ergodic quotients (bottom) are given 
from left to right for a harmonic oscillator, a Dufhng-type 
oscillator and a double gyre system. Dots indicate equi- 
libria and dashed arrows indicate how moving an initial 
condition in the state space moves its image on the er- 
godic quotient. 

In Figure [5] we demonstrate what the ergodic quotient 
might look like locally for several simple state-space por- 
traits. In elliptic regions, invariant functions are a one- 
parameter family, mapping to a filament in the ergodic 
quotient. In DufHng-type dynamics, one-parameter invari- 
ants exist along libration trajectories in elliptic zones and 
along the revolution trajectories outside of the separatrix. 
Moving initial conditions close to separatrix from either 
side results in trajectories slowing down more and more as 
they pass the saddle, and consequently time averages of an 
observable will converge to the value that the observable 
attains at the saddle. As a result, filaments for librations 
and revolution trajectories meet in the ergodic quotient. 
In double-gyre flows, the filaments do not meet, since two 
adjacent elliptic cells share only two out of four saddles; as 
trajectories slow down next to all four, the time averages 
will converge to a weighted average of the function values 
at the saddles, which will be different for every adjoining 
elliptic cell. 



Remark. Fomcnko in [2[ used a similar construction, 
the Reeb graphs, to analyze integrable systems. Reeb 
graphs are graphical representations of connected compo- 
nents in level sets of Morse functions, in particular, the in- 
tegrals of motion. For an integrable system, the joint level 
sets of integrals of motion are ergodic sets. Fomenko did 
not extend his work to systems for which we do not have 
expressions for integrals of motion. The ergodic quotient, 
therefore, can be interpreted as a setting for generalization 
of Fomcnko's approach to non-integrablc systems. 

2. 4- Diffusion modes as global coordinates 

Although the ergodic quotient is a subset of an infinite- 
dimensional sequence space, its main features could be 
approximated well in a lower number of dimensions, as 
illustrated in Figure [5] To obtain a low-dimensional rep- 
resentation, we will convert from the ambient coordinates, 
where axes are given by averaged observables, to intrin- 
sic coordinates, where axes are independent parameters 
that parametrize the ergodic quotient. As intrinsic coor- 
dinates, we chose the modes of diffusion, or heat spread, 
along the ergodic quotient, which provide an efficient, or- 
thogonal low-dimensional representation. Numerically, we 
pass from ambient coordinates to diffusion modes using 
the Diffusion Maps algorithm Q. 

Problems of approximating high-dimensional data by 
low-dimensional parametrizations are common in data min- 
ing and machine learning communities, where a low-di- 
mensional process, e.g., an object moving in the field of 
vision, is measured by high-dimensional observation, e.g., 
image consisting of pixels. The usefulness of diffusion 
eigenfunctions, as a prominent example of laplacian-based 
methods, is well documented in such settings, e.g., j4fj| . 

On a differentiable manifold, the diffusion operator is 
given by A = e~ tA , where A is the Laplace-Beltrami oper- 
ator and t the time over which diffusion evolves. Spectrum 
of the diffusion operator is contained within the unit in- 
terval, with first eigenvalue always Ao = 1, and the rest 
of the eigenvalues ordered in decreasing order. Eigenfunc- 
tions of the diffusion operator x^ are orthogonal and nor- 
malized to ||x amplitude with the trivial eigenfunc- 
tion = 1- The numerical algorithm, The Diffusion 
Maps, retrieves the values of eigenfunctions sampled at 
finite number of points on the underlying manifold. 

When we apply Diffusion Maps to sequences of time 
averages, the ergodic quotient plays the role of the differ- 
entiable manifold. The interpretation of numerical results 
as discretizations of the Laplace-Beltrami eigenfunctions 
is strictly valid only in the case when ergodic quotient is a 
differentiable manifold, as it seems to be for some simple 
systems, but not necessarily in the general case. Despite 
this fact, this paper demonstrates that Diffusion Maps are 
practically useful in analysis of dynamical systems and we 
expect that a consistent interpretation will be found even 
for sets without a differential structure. 

The properties of diffusion eigenfunctions are demon- 
strated on simple examples in papers 0, , and detailed 
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analysis of the mathematical properties can be found in 
18, 2(J 38 1 . In short, we can expect the lowest-index eigen- 
functions to be combinations of indicators on disconnected 
components of the ergodic quotient. The eigenfunctions of 
higher index behave similarly to Fourier modes on each of 
the components: first varying over coarser features and, as 
the order is further increased, varying over finer features. 

In the spirit of representation sequence © , we will add 
embedding into the eigenfunctions of the diffusion operator 
as the final step: 

X A fJ, x ^ (. . . , f k (x), . . . ) k£Z D 



in 



(X lX {1 \x),X 2X ^(x),...), 



and refer to embedding x i— > (XiX ( x )> \2 X ^(x), • ■ • ) as 
the diffusion coordinates of the dynamics. The sketch of 
the embeddings was given at the beginning as Figure [TJ 



diffusion distance 




Figure 3: A sketch of distances in the space of trajectory 
averages. The full line represents the ergodic quotient, the 
chordal distance is the empirical metric, and the dotted 
line indicates the point where a small amount of noise can 
change the geodesic distance discontinuously, but not the 
diffusion distance. 

The cuclidcan distance in the diffusion coordinate rep- 
resentation carries a geometric meaning, which is different 
from, but related to the meaning of the distance in the 
ambient space of the ergodic quotient, spanned by aver- 
aged observables. The empirical distance T> is the chordal 
distance in the ambient space, i.e., the length of a short- 
est path between points. Diffusion distance, like geodesic 
distance, takes into account only paths that lie within the 
ergodic quotient. However, while the geodesic distance is 
the length of the shortest path connecting two points, dif- 
fusion distance is, in a sense, the average length of all paths 
that connect the two points. This average is computed as 
the difference in heat flow reaching two points when heat is 
injected at any other point on the quotient, averaged over 
all possible positions of the heat source. A more detailed 
description can be found in 0] . 

For our purposes it is sufficient to state that averaging 
over all paths provides robustness to noise. If the noise 
introduces a short-cut between two parts of the manifold, 
as depicted in Figure [31 the geodesic distance changes dis- 
continuously, however, the diffusion distance does not, as 
averaging over all paths ensures that a small amount of 
noise does not change the geometry significantly [13; §4-2]. 

The Diffusion Maps algorithm approximates the op- 
erator e~ tA by a N x N matrix, where N is the num- 



ber of points sampled from the ergodic quotient: in our 
case equal to the number of initial conditions. Discrete 
diffusion modes are obtained as eigenvectors of that ma- 
trix. The |Appendix C| provides more details about the 
implementation: the input to Diffusion Maps is a matrix 
of pairwise P_ s distances between computed trajectories, 
while the output is the set of diffusion coordinate vectors, 
assigning diffusion coordinates to each trajectory. 

2.5. Example: Steady ABC Flow 

We will demonstrate the use of the ergodic quotient on 
the Arnold-Beltrami-Childress (ABC) flow 0. This flow 
was previously used to demonstrate identification of La- 
grangian Coherent Structures in (l4l ] and almost-invariant 
sets in [l2| ■ It is a steady, volume-preserving flow, evolving 
on a 3-torus. Despite its very simple system of ODEs, the 
state space of the ABC flow contains both regular vortices 
and chaotic zones. 

Trajectories of the ABC flow obey the following equa- 
tions on (x, y, z) e T 3 



A(x,y,z) = y 



A sin 2irz + C cos 2ny 
B sin 2irx + A cos 2ttz 
C sin 2iry + B cos 2irx 



(5) 



As in [IJ, [lj], we used A = V3, B = V2, C = 1, but 
the domain was rescaled to VolT 3 = 1. Observables are 
the complex Fourier harmonics (UJ with wavevectors k € 
[— 10, 10] 3 n Z 3 . The trajectories were integrated using 
an adaptive solver (Radau) for at least T = 500 per ini- 
tial condition, but possibly longer to achieve tolerance 
ATOL = 2 x 10~ 4 in convergence of averages for each 
observable and trajectory (see | Appendix A| q 

In short, the steps used for identification of coherent 
structures arc: 



1. Select initial conditions 



2. Compute trajectories and time averages until con- 
vergence, 

3. Compute H~* distance matrix 

(dij) = V- s (fi Xi ,fi X] ), 

4. Feed (dij) into Diffusion Maps to compute diffusion 
coordinates (Xi X i(xi), ^2X2(%i), ■■■), 

5. Compute k- means clusters (optional, for visualiza- 
tion) . 



3 The precise value of the tolerance chosen did not make a big 
difference in data obtained by the method for the ABC flow. Reduc- 
ing tolerance further extends the computation on a relatively small 
number of trajectories (see Figurc lA.911 while most of the trajectories 
require very similar simulation times even if the tolerance is reduced 
by an order of magnitude. 
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(a) Six subsets in the state space corresponding to six clusters 
of aggregation. Subsets contain primary vortices of the ABC 
flow. The chaotic sea between vortices is the seventh cluster 
and appears as the void between vortices. 



Index 4 -°- 08 -0.15 Index 2 

(b) Projection of ergodic quotient onto three out of 
ten diffusion coordinates used for clustering, axis la- 
bels are indices k of coordinates A^X ■ Points were 
colored according to membership in clusters. Cluster 
at the top of the figure corresponds to the chaotic sea. 



Figure 4: Six primary vortices extracted by k-means clustering (k = 7) of projection of ergodic quotient onto first 10 
diffusion coordinate. ABC flow ([5]) was simulated with A = B = \/2, C = 1, from N = 1000 initial conditions 
uniformly distributed in basic periodicity cell [0, l) 3 , with observables cut off at wavenumber K = 10, and convergence 
tolerance ATOL = 2x 10" 4 . 



Connected components in the ergodic quotient corre- 
spond to sharply distinct coherent features. Information 
about all distinct components will rarely be visible from 
a single, or a few, diffusion coordinates. To quickly re- 
cover the dominant features in the flow, we can apply a 
clustering algorithm. Clustering based on the euclidean 
distance, e.g., the fc-means algorithm, is justified since the 
euclidean distance approximates the diffusion distance on 
the ergodic quotient. 

The Figure 0] demonstrates the application of k-means 
clustering, with k = 7, to the first ten diffusion coordi- 
nates on the ergodic quotient of the ABC flow. The six 
clusters correspond to six primary vortices that were iden- 
tified in 0, Fig. 4]. Furthermore, these are the same 
structures visible in LCS analysis of the flow 14j, Fig. 1] 



and the eigenvectors of the Ulam discretization of Perron- 
Frobenius operator [l2|, Fig. 6]. This corroborates our 
assertion that the components of the ergodic quotient cor- 
respond to features typically identified as coherent struc- 
tures in fluid flows. 

In contrast to LCS and Ulam's methods, the initial 
conditions do not have to be placed in the entire state 
space, as mentioned in Section [2. II it is sufficient to place 
them on a set that intersects interesting coherent struc- 
tures that we want to analyze or visualize. In Figure [S] we 
show a visualization based on N = 500 initial conditions 
sourced from a rectangle (y, z) £ [0.35,0.8] x [0.6,0.9], on 
the x = face of the unit cell. This region intersects one 
of the six primary vortices shown in Figure I4al Since we 
focused the initial conditions in a smaller region of state 
space, the same number of diffusion coordinates discerns 



smaller features in the state space, in this case the sec- 
ondary 1 : 2 vortices that wrap around the central vortex. 

The Figure [5c] shows that some diffusion coordinates 
are local-ordering functions on the elliptic zones, acting 
as computationally constructed local integrals of motion. 
The 10th ergodic coordinate acts as an energy-like inte- 
gral of motion, locally on the inner primary vortex. The 
blue filament in upper-left corner of Figure I5bl indicates 
that the positive values of the 12th coordinate vary ap- 
proximately continuously over the secondary vortex lobes, 
colored in blue in Figure I5al This is an indication that a 
plot similar to Figure [5c] could be made to demonstrate 
local-ordering on a part of the secondary vortex, using the 
12th coordinate. 



3. Extensions including time- varying features 

3.1. Periodically-driven flows 

Although we have presented the construction of the er- 
godic quotient for the autonomous flows, in practice, noth- 
ing prevents one from applying the exact same computa- 
tional procedure to a non-autonomous flow. In this Sec- 
tion, we explain that when system is non-autonomous, yet 
periodically-driven, one still obtains useful and non-trivial 
information about the structure of invariant sets. 

The non-autonomous systems can be embedded in an 
autonomous system by extending the state space. Let 
Ai e = Ai x T be an extension of the state space Ai. 
For periodically driven flows, T can be taken as T, while 
for aperiodically driven flows and finite-time flows T can 
be taken as R and [ti,tf] C R, respectively. 




(a) Five subsets identified by the k-means algo- (b) Projection of the ergodic quotient (c) Level sets of the 10th diffusion eigenvector, 
rithm with k = 5 acting on first ten diffusion onto three diffusion coordinates, axis la- The chaotic sea has been cropped for clarity of 
coordinates. bels are indices k of coordinates \kX^ ■ display by thresholding the first diffusion eigen- 

Points were colored according to mem- vector. 

bership in clusters. 

Figure 5: A primary vortex aligned with the x-axis of the ABC flow (A = \^3, B = y/2, C = 1). Lobe-shaped structures 
that wrap around the central vortex correspond to secondary 1 : 2 vortices. Results are based on N = 500 trajectories 
initialized uniformly in rectangle x = 0, (y, z) G [0.35, 0.8] x [0.6, 0.9] , with obscrvables cut off at wavenumber K = 10, 
and convergence tolerance ATOL = 2 x 10~ 4 . 



For a non-autonomous dynamical system defined by 
an ODE x = A(x,t), where i e M, we can define the 
extended, autonomous system through 

x = A(x, t), 
f = c, 

where [x, r) S A4 e and c is a time-rescaling constant. For 
periodically-driven systems on compact M., i.e., for T = T, 
the extended state space M e is also compact. If, addition- 
ally, the flow map ip t : M e — > M. e of the system is continu- 
ous, the system (M e , (fit ) fulfills all the conditions that we 
required for the construction of the ergodic quotient map 
7r e : M e -> given in Section |2~T1 

The difference between an autonomous system on M. 
and a periodically-driven system cast into an autonomous 
system on M e is in the eye of observer. Typically, we 
will treat only the first factor of M e as "physical" states, 
always keeping in mind that T was a formal construct 
accounting for the "physical" time. This is reflected in 
observables usually still available as functions / : M. — > C. 

If the observables on M. were chosen as Fourier func- 
tions fk : M. — > C, with k £ Z D , we can interpret them as 
a subset of Fourier observables g w : M e — > C, w € Z d+1 , 

g w (x,r) =Cf k (xy 2 ™^ T , 

with the constant C normalizing the observable to Hfl^H^ = 
1. When wave-vector w is constrained to w = (fc,0), the 
identity g w {x,T) = fk{x) holds for any r. Consequently, 
the trajectory averages g w (x,r) = fk(x) are then insensi- 
tive to the initial "physical" time r. 

If we proceed with such a partial set of observables g w , 
where w = (k, 0), k G IP ', we cannot in general hope to 
obtain the full ergodic quotient by mapping x — > g w (x, r). 



Denote the partial quotient map by tt\m ■ M e — > l°°(Z D ), 
defined through 

k\m{x,t) :=(... ,g w (x,r), ... ) w , 

for w = (k,0) and k £ Z , with the associated partial 
quotient £\m = t^\m Looking at pre-images of any 
trajectory-averaged observable, it follows that 3^ 1 (c) = 
/ fc _1 (c) x T. Therefore, for any point p € IP , with pre- 
image 7r _1 (p) formed using observables fk, the pre-image 
under partial quotient map is tt\m(p) = tt^ 1 (p) x T. Such 
pre-images are still invariant sets for the dynamics, as they 
are formed from joint level sets of averaged observables. 
The difference between the invariant partition correspond- 
ing to the ergodic quotient of the extended system £ e and 
the invariant partition corresponding to the partial quo- 
tient is that only £ e yields the ergodic partition, in 
general. 

In plainer terms: any analysis based on the partial quo- 
tient £\m will identify features distinct only in physical 
state-space M , but not distinctions in physical time. This 
is sometimes desirable; as a trivial example consider a sim- 
ple oscillator, 6 = 1 for 6 € A4 = T, and, even though it 
is already autonomous, consider its formal extension into 
A4 e = T 2 . Analysis of the partial quotient £\m would iden- 
tify a single invariant set M. , while the analysis of £ e would 
identify a continuum of ergodic sets, each corresponding 
to a trajectory with a different initial phase. Depending 
on the application, either of the two analysis can contain 
valuable information about the analyzed system. 

A common technique used in analysis of periodically- 
driven systems is the Poincare map $™ :X->jVl associ- 
ated to the system: 
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It is a simple result to show that for continuous flows 
ip t , the choice of the initial time (or phase) to results in 
topologically-conjugate maps, implying that each there ex- 
ists a homotopy such that each invariant set of $ tl has a 
homotopy-related invariant counterpart in state space of 
$t 2 . However, the analysis presented in this paper goes 
beyond topology, as we study geometry of the features in 
the state space. 



,4 



B 



o 



2tt 



Figure 6: Three trajectories illustrating the distinction be- 
tween the direct study of a Poincare map and coherent 
features identified by partial ergodic quotient analysis. 

Consider three trajectories sketched in the Figure H)] If 
a Poincare map was taken at t = (mod 2ir) all of them 
would be represented as fixed points in state space of the 
map. Looking at just the Poincare map, we would be much 
more likely to identify trajectories A and B as similar, 
since they pierce the Poincare surface at nearby locations. 
However, considering the entire evolution, trajectories B 
and C deviate significantly only over relatively short time 
interval. Since time averages capture mean behavior, H~ s 
distance on tt\m would (correctly) identify B and C as 
similar, regardless of the time instance at which Poincare 
section used for visualization is taken. 

In summary, constructing the partial quotient, i.e., em- 
bedding the state space using trajectory averages of ob- 
servablcs on the "physical" states Ai, identifies ergodic 
partition of the Poincare map. Additionally, the met- 
ric structure on the partial quotient reveals which invari- 
ant sets are similar based on the entire trajectory evolu- 
tion, not just on the points at which trajectories pierce 
the Poincare surface. Therefore, this method extends the 
Poincare map analysis by providing additional information 
which is potentially useful in identifying coherent struc- 
tures. To demonstrate, we proceed with a numerical ex- 
ample. 

3.2. Example: Periodically- driven 3D Hill's Vortex Flow 

As an example of a periodically-driven flow, consider 
the following ODE, evolving states x = (R,z,9) C M + x 
IxT: 



R 




2Rz 


z 




1-4R- 


6_ 




c/2R 



V2i?sin6» 
z(V2R)- 1 smt 



2 cos 6» 



sin27ri, 



where parameters c and e are swirl and perturbation strengths, 
respectively. For all values of c and e the system pre- 
serves the state space volume, i.e., V • A = 0, where 
V = (dR,d z ,R~ 1 dg). This system was studied analyti- 
cally in 29|, [39| as an example of a non-hamiltonian system 
showing Kolmogorov- Arnold- Moser-type structures at low 
values of e = c. In this paper, we will show the struc- 
ture of invariant sets for both low- and high-valued per- 
turbations for the purposes of illustrating the method of 
analysis; therefore, we will not explore here the bifurca- 
tion mechanism behind the structures in the state space. 
A preliminary numerical study was presented in Q ■ 

The results of the first simulation (Figure [7|) demon- 
strate the low-perturbation structure of the state space, 
at c = e = 0.01. Uniformly distributed initial conditions 
(N = 1000) were seeded on the 9 = plane in the rectangle 
(R, z) e [0.01, 0.30] x [-0.5, 0.5], which is within a bounded 
invariant set. The observables were chosen as Fourier ba- 
sis functions on [0.0,0.5] x [-1.0,1.0] x [0,27r] C M, with 
wavevectors in the box \kn\ < 8, \k z \ < 8, \kg\ < 4. Sim- 
ulations were run for no less than T m j„ = 600, and until 

was achieved (see |Appendix A 
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A(x,t) 



tolerance ATOL 

for definition), with most simulations terminating before 
2000 time units expired. 

At c = e = 0, any fixed-0 plane is invariant and on it 
the system conserves the Hamiltonian H (R, z) = i?z 2 — R+ 
2R 2 . The Figure [7] demonstrates that at low perturbation 
values, the structure of the cross-section of invariant sets 
is very similar to the level sets of the Hamiltonian, despite 
the flow having a full three-dimensional character. At the 
plotted scale, \i appears to be a monotonic function of the 
Hamiltonian; however, the non-hamiltonian KAM theory 
predicts that invariant tori do not form a continuum, but 
are interspersed with thin chaotic zones. Observables of 
low wavenumbers effectively coarse-grain dynamics; since 
the chaotic zones are so thin at low perturbation values, 
they are not identified as significantly different invariant 
sets. Plotting higher-order diffusion eigenfunctions against 
the xi reveals the similarity of diffusion eigenfunctions 
to the Legendre polynomial system, which is likely due 
to connection of Legendre polynomials to spherical har- 
monics and the spherical shape of the Hamiltonian func- 
tion. In this case it is particularly obvious how mono- 
tonic parametrization of the ergodic quotient can serve 
as a proxy for conserved quantities, possibly leading to 
Fomenko-type analysis. 

When the perturbation is increased, the dissolution of 
invariant tori continues, but an additional invariant set be- 
comes visible in the state space. For the next simulation, 
we have chosen c = e = 0.3495. There is nothing special 
about this particular parameter value: all perturbations 
c = e between roughly 0.2 and 0.4 exhibit the same topo- 
logical structure, visible on similar scales. As the region 
where trajectories stay bounded shrinks with the growth 
of the perturbation, we have seeded the initial conditions 
uniformly (N = 1000) in (R, z) € [0.2,0.3] x [-0.1,0.1] in 
9 = plane, with initial time t = 0. The observables were 
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(a) The slice at 9 = through the Poincare section at (b) Points of the Poincare section used for interpo- 

t = (mod 1) in the invariant region; color interpolated lation in panel [(a)] colored using the first diffusion 

using the first diffusion eigcnfunction xi from points in eigenfunction xi- 
panel |(b)| 




0.125 0.25 0.375 0.5 

R 



(c) Bounded level sets of the unper- 
turbed Hamiltonian H(R, z) = Rz 2 — 
R + 2R 2 in a (R, z)-plane. 
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Xi 

(d) Diffusion eigenfunctions Xfc, k = 
2, 3, 4, 5 of the ergodic quotient plot- 
ted against xi- 



Figure 7: Structure of the invariant sets in the state space of a periodically-forced 3D Hill's vortex at low perturbation 
c = e = 0.01. 



Fourier functions on [0.0, 0.5] x [—0.5, 0.5] x [0, 2ir), with the 
same wave-vector bounds as before. The simulation times 
and tolerances also remained as in the low-perturbation 
case. 

The Figures [5a] and [8b] show a secondary coherent set 
that becomes visible within the onion layers at higher 
perturbation values. The primary tubular invariant set 
(in red, Cluster 3) remains at the center, with the sec- 
ondary invariant set (in blue, Cluster 0), sickle shaped in 
its cross section, wrapping around the tubular core in 1:2 
resonance. Both primary and secondary sets are enclosed 
by the outer tubular shells (in yellow, Cluster 2). Here, 
we stress again, that even though we are showing sets in 
a three-dimensional Poincare section, they are truly time 
invariant, despite the time-dependent forcing. More pre- 
cisely, to each of the invariant sets C C M identified by 
clustering, there exists a direct extension in the extended 
state space CxTcM e . 

The topology of the embedded quotient reflects the 
change in the invariant set structure, compared to the low 



perturbation case. While embedding in xi vs - Xi i n Figure 
I7dl showed a topological line interval, the high perturba- 
tion case in Figure [8c] shows that a branch splits off the 
interval which was once formed by segments correspond- 
ing to clusters that compose the inner and outer tubular 
regions. The Figure [8d] takes a closer look at diffusion 
eigenfunctions within only one of the identified clusters. 
Comparing the shapes of coordinate functions with the 
low perturbation case in Figure 1751 in interval \i <S [0-4, 1], 
we can again notice the similarity to Legendre polynomi- 
als, indicating that the flow appears to locally conserve 
a continuous motion integral, even if that is not the case 
globally. 

3.3. Harmonic quotients 

While this paper deals primarily with time-invariant 
structures, the Section ^. ll describcs how even time- varying 
flows can be analyzed using trajectory averages, through 
embedding of the system into an extended state space, 
expanded by time variable. Fourier observables g on the 



11 




0.1 



(a) The Poincare section at t = (mod 1). 
Segments of state space corresponding to clus- 
ters and 3 shown in blue and red, respec- 
tively. The remaining clusters are cropped for 
clarity. 
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(c) Embedding of the partial quotient 
into first two diffusion eigenfunctions 
Xi> X2; with color indicating the clus- 
ter assigned to the point. 
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(b) The slice at 9 = through the Poincare sec- 
tion in panel |(a)| displaying all clusters; colors 
reflect cluster labels (0, 1, 2, 3). 




(d) Diffusion eigenfunctions \k °f the partial quo- 
tient, plotted against \l f° r points in the cluster 
2 



Figure 8: Structure of the invariant sets of a periodically-forced 3D Hill's vortex at high perturbation c = e = 0.3495. 
Diffusion coordinates were clustered using fc-means algorithm into 4 clusters to extract the main features. The Poincare 
section used for visualization is at t = (mod 1). 



extended state space M e = M x T were formed from 
Fourier observablcs f k : M — s- C as g w (x,t) = fk(x)e l2nut , 
where w = (k, w) was the extended wave- vector, with u> the 
wavenumber along time axis. Time averages of observables 
were given by 

1 f T 

9w(x,t)= hm — / g[ip T (x),T]dT 

= Hm I f[fo VT ]{x)e a —dr. 

T->oo T J 

When T = T, it was sufficient to consider only uj G Z 
to obtain a full function basis on M. e . However, when the 
time is the full real line, as in the aperiodic case, we require 
a real- valued wavenumber The complete functional 

basis onM e is then uncountable, and the crgodic quotient 
of the extended system, evolving in A^ e , would not be a 
subset of a sequence space anymore. 

Nevertheless, even in the case of flows cpt generated 



by autonomous systems, it is of interest to fix lu and to 
consider functions given by averages 

f^\x) 4 lim i f T [f o <p r ]{ X y % ™ T dT, 

termed harmonic averages. The harmonic averages were 
considered by Wiener and Wintncr under the name Fourier 
averages [4l|, but more recently, they have been analyzed 
in [27[. We use them here to provide a generalization of 
the ergodic quotients which can be used to analyze peri- 
odic and quasi-periodic transport. 

For any fixed u) € R, and Fourier observables fk, the 
embedding tt u (x) — (. . . , fj^'' \x), . . . )k£z D generates an 
analog of the ergodic quotient, the harmonic quotient at 
frequency lu. The ergodic quotient is then a special case 
of the harmonic quotient, computed for co — 0. In under- 
standing the crgodic quotient, the level sets of averaged 
observables / played a crucial role as invariant sets. Sim- 
ilarly, level sets of provide insight into periodic sets 
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and periodic transport when ui € Q, and wandering sets 
and quasi-periodic transport when u g Q [27j|. More con- 
cretely, if w = 1/2, will vanish over all points x that 
are not (^-periodic with period 2. Consequently, the joint 
level sets of away from (0, 0, ... ) identify sets in state 
space between which there exists dynamical transport of 
period 2. 

The functional mapping / i— > f( u ' has an interpretation 
connected to the Koopman operator. For systems satisfy- 
ing conditions in this paper, the Koopman group consists 
of linear, bounded composition operators Ut ■ T — > J-, 
defined as 

Utf - f°Vt, 

for time t € T. The space J 7 is a space of observablcs, 
in literature typically chosen as L°° (Ai , fx) , where fj, is a 
^(-invariant measure. The Koopman operator is the gen- 
erator of the Lie group; when time is discrete, i.e., TcZ, 
the generator is often denoted simply by U. 

The harmonic averages f^> are eigenfunctions of the 
Koopman operator, since it is simple to show that 

Utf^ — e i27TUJt f( u \ 

It follows that, for any observable /, the harmonic aver- 
age is an eigenfunction of the Koopman group gen- 
erator at eigenvalue A w = e l27r " for discrete time, and at 
eigenvalue A w = i2nu} for continuous time. The operator 
P w : T -> T, 

P u f = f (u) 

is a projection operator, projecting the space of observ- 
ables onto the eigenspace of the Koopman operator at A w . 
In this context, the harmonic quotient map 7r u , applied 
to the Fourier basis, evaluates the functions spanning the 
eigenspace at A w . Therefore, by endowing the harmonic 
quotient with a particular geometry, we are able to compu- 
tationally analyze the cigenspaces of the Koopman opera- 
tor, without ever computing, or approximating, the Koop- 
man operator itself. 

4. Discussion 

It is difficult to give a level-field comparison of different 
methods for identification of "coherent structures" in state 
spaces due to the lack of a method-independent, universal 
definition of what a coherent structure might be. While 
each of the methods, namely, LCS, Ulam's approximation, 
and ergodic quotient, has an operational understanding of 
a coherent structure, ultimately the experience and appli- 
cation is the only judge of how well a method captures 
what is identified as "coherent" in practical settings. For 
this reason, we decided to apply the analysis to the ABC 
flow; this was not because we viewed the ABC flow as 
a particularly challenging problem, but rather because it 
was well analyzed both theoretically and computationally, 



resulting in well-established knowledge about what coher- 
ent structures one should obtain. In this sense, we found 
that the analysis of ergodic quotient identified the same 
structures as the theoretical analysis and the alternative 
computational methods. 

When it comes to the amount of work, computational 
or experimental, needed to compute an approximation to 
invariant sets, one has to asses the setting in which either 
of the methods is applied. An advantage of the ergodic 
quotient method is that it requires a single initial condi- 
tion to be placed in any ergodic set sought to be identi- 
fied. In practical terms, this means that simulations or 
experiments need to be initialized only on a small sub- 
set of the entire state space: this was clearly exemplified 
in Section 13.21 where simulations were initialized on the 
9 = planes and the method successfully identified coher- 
ent structures in the full (R, z, 9) state space. By contrast, 
methods based on Ulam's approximation require seeding 
the initial conditions in the entire invariant domain. While 
this might not be a significant distinction when working 
with ODEs, in experimental settings, e.g., when analyzing 
fluid flows using Particle Tracking Vclocimetry, it might 
be possible to initialize the trajectories only in a specific, 
small region, and count on the long experimental runs to 
collect data about the rest of the state space. Such data 
is a natural input to ergodic quotient methods, whereas 
Ulam-type computations would require additional process- 
ing, e.g., flow interpolation, that might introduce further 
errors. 

The number of initial conditions required to approxi- 
mate the ergodic quotient is fairly low: only TV = 1000 
trajectories were used to retrieve a very detailed image 
(S = 200 3 spatial cells) of the flow structures in Section 
12.51 This is a consequence of the recognition that the entire 
trajectory can be described to a high precision by a rela- 
tively small set of values, i.e., averages of a function basis. 
While the entire trajectory evolution is used for visual- 
ization of the structures, the identification of structures, 
where the algorithmic effort is really spent, is performed 
only with TV data-points (samples of the ergodic quotient) 
as inputs. 

In terms of numerical linear algebra, we do acknowl- 
edge that although the Ulam's matrix is large, S x S, 
where S is the total number of spatial cells, it is often 
sparse. In contrast, the eigenproblem at the core of the 
Diffusion Maps is is potentially dense, although on a much 
smaller TV x TV matrix, where TV is the number of trajecto- 
ries]! Nevertheless, we see the distinction as conceptual, 
rather than computational, signifying that the complex- 
ity of the state space is often captured by a far smaller 
number of parameters than needed to describe transport 
between pairs of cells in state space which is especially true 



4 Due to the rapidly-decaying exponential diffusion kernel, how- 
ever, the diffusion eigenproblem might also be sparsified by truncat- 
ing the kernel. We have not explored the consequences of sparsifica- 
tion, though. 
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for (near-) integrable, high-dimensional systems. This pa- 
per focuses on visualization as the first step in validation 
of the ergodic quotient approach; however, we expect that 
the ergodic quotient could be analyzed directly, for exam- 
ple, to detect topological changes such as those between 
embedding \i vs - X2 in Figures ITdl and l8cl These changes 
indicate a bifurcation of the flow, without ever visualiz- 
ing the state space. If our expectations are confirmed, the 
high-dimensional state spaces, which cannot be visualized 
and requiring a large number of cells to discretize the state 
space, could potentially be analyzed through the study of a 
small number of trajectories by noting topological changes 
in the ergodic quotient. 

Both the Ulam's approximation and the ergodic quo- 
tient construction feature exponential growth of problem 
size with the increase of the state space dimension D. 
Number of state space discretization cells, for a fixed axial 
resolution R, grows as R D . Similarly, to approximate the 
ergodic quotient, if we truncate the function basis at a uni- 
form cutoff wavenumber K, the number of required Fourier 
obscrvables is K D . For certain high dimensional systems, 
it is not always necessary to analyze the state space di- 
rectly: rather a particular, low-dimensional output space 
might suffice. For example, in networks of power genera- 
tors the variable of interest in detecting instabilities might 
be the mean frequency across all generators (see [36|). If 
instead of the state space we study observables only on the 
output space, the number of required observables can be 
significantly lower. It is therefore possible that it is suffi- 
cient to initialize the trajectories on a low-dimensional set 
in the state space and analyze it in a low-dimensional out- 
put space, despite the dynamics being high-dimensional. 
Both order reductions are a natural fit with the presented 
technique in analysis of trajectory averages presented here. 

The application of the ergodic quotient techniques to 
periodically-forced systems has been presented in Section 
[21 The extension of the method to the case of transient sys- 
tems, where long-time averages lose the information about 
the initial evolution, is the focus of our future studies. The 
case of transient systems is related to the study of finite- 
time systems, i.e., those defined only on a fixed time in- 
terval, which includes all data sets where experiment time 
cannot be extended into the steady regime, e.g., field sam- 
ples of geostrophic flows. We do not intend delve deep 
into a discussion about how much we can predict about 
the future of an aperiodic flow based only on a finite time 
interval. Instead, we point out that the analysis presented 
in this paper does have a very clear interpretation even on 
a fixed finite-time interval: the finite-time quotient would 
reflect sets of initial conditions that remain together under 
advection by the flow map over the interval considered. In 
that case, the H~ s metric and diffusion maps result in a 
multi-scale aggregation of finite-time trajectories, identi- 
fying regions of space where material is transported in a 
homogeneous fashion. The resulting "coherent structures" 
present a template of the transport of material over the an- 
alyzed time interval. It is left for future research to study 



the relationship between transport templates correspond- 
ing to overlapping or adjoining time intervals and their 
invariance and predictive power. 

Recently, the analysis of mesohyperbolicity and mesoel- 
lipticity has yielded success in the study of finite-time 
flows. Both the ergodic quotient and mesohyperbolicity 
depend on averages of observables along trajectory lines. 
The first distinction between these methods is in their pur- 
pose: mesohyperbolicity and mesocllipticity arc proposed 
generalizations of concepts of hypcrbolicity and cllipticity 
from autonomous, steady-state flows to unsteady, finite- 
time flows. In contrast, ergodic quotient is directly related 
to the ergodic partition, which studies invariant sets. The 
two concepts are potentially related, as is the case in clas- 
sical theory of autonomous systems, where the relation- 
ship between analysis of invariant manifolds and analysis 
of Lyapunov coefficients is standard. 

The second, more technical, distinction is in the objects 
averaged along trajectories. In construction of the ergodic 
quotient, the ultimate goal is to include all observables 
in the averaging process, in an attempt to represent the 
underlying empirical measures as well as it is possible. The 
observables chosen depend only on the topology of the 
state space, not on the actual flow map of the dynamical 
system. In studies of mesohyperbolicity /ellipticity, a very 
particular set of observables is chosen: they are entries 
in the Jacobian matrix of the flow. Therefore, while the 
methods use the same technique of computing averages of 
obscrvables, the theory behind the averages serves different 
purposes, and, ultimately, demonstrates the versatility of 
information that can be extracted from trajectory averages 
of observables. 

From the standpoint of the direct application of meth- 
ods presented in this paper, there arc still points where 
one could improve the presented procedure, although most 
of such efforts fall outside the scope of dynamical sys- 
tems, e.g., into fields of machine learning or computa- 
tional geometry For example, setting the bandwidth of 
the exponential kernel in the Diffusion Maps algorithm 
(sec I Appendix Cj ), although relatively simple to perform 
manually, is still not fully automatic, which might slow 
down the analysis process when dealing with a particu- 
larly complex ergodic quotient. 

Furthermore, we were guided mostly by simplicity and 
clarity in our choice of the fc-means algorithm for post- 
processing of diffusion coordinate, in an effort not to de- 
rail the main point of this paper. The Figure I5cl reflected 
the imperfections in the choice of fc-means clustering as 
the final step in post-processing. One could expect that 
the four clusters identified would be the three "branches" 
and the point at which they merge. While the numerical 
results come close to such classification of points, there is 
room for improvement. It is possible that instead of iden- 
tifying clusters, identification of connected components or 
one-dimensional sets, based on a proximity graph in a dif- 
fusion coordinate space, would yield even better results 
than those presented here. 
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5. Conclusion 

Our results demonstrate that, by analyzing the ergodic 
quotient, we can infer how trajectories organize into coher- 
ent structures in the state space. To construct an approx- 
imation to the ergodic quotient, we average a subset of a 
function basis on the state space, and endow it with the 
empirical distance. In this setting, the geometry of the 
ergodic quotient reveals dynamically coherent regions in 
the state space. The geometry was studied using diffusion 
coordinates which rectify the ergodic quotient, enabling 
clear visualization and extraction of clusters correspond- 
ing to coherent structures. 

To evaluate the method, we identified coherent struc- 
tures for the ABC flow, an area-preserving, three-dimensi- 
onal, steady flow. Our results correspond to those coherent 
structures identified by Lagrangian Coherent Structures 



14| and Ulam's method (I2j. Unlike those approaches, 



however, our method does not require placing the initial 
conditions in the entire state space, but rather depends 
on long-time integration to explore the state space. For 
the same reason, it requires a relatively small number of 
initial conditions. This makes it suitable for analyzing ex- 
perimentally harvested trajectories in experiments where 
initial conditions can be placed in a small region of the 
state space, but the dynamics of the system evolves on a 
larger domain. 

For autonomous flows, points on the ergodic quotient 
represent ergodic invariant sets in the state space. When 
the flow is periodically-dependent on time, the construc- 
tion we present identifies invariant sets of the Poincare 
map. In contrast to analysis of Poincare map, however, 
coarse graining of invariant sets takes into the account the 
entire evolution of the trajectory, not only their Poincare 
section, providing a more accurate method of identifying 
trajectories which have similar behavior. On an exam- 
ple of a periodically-forced 3D Hill's vortex flow we show 
that the method is able to identify previously undescribed 
coherent features in state spaces. Finally, we close the 
discussion of time-dependent features by generalizing the 
ergodic quotient to harmonic quotients, which can provide 
information about periodic and quasi-periodic transport 
in the state space. 

Although this paper discusses applications of the er- 
godic quotient to analysis, we believe that the proposed 
framework could, in future, be useful in design. Diffu- 
sion coordinates provide an orthogonal system of invariant 
functions for the system. They are, therefore, particularly 
suitable as a setting for design of invariants of motion, 
e.g., through optimization of linear combinations of diffu- 
sion coordinates based on a design criterion. 
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Appendix A. Convergence of averages 

Our implementation of the averaging process computes 
averages of obscrvablcs along with the numerical evolu- 
tion of the trajectory, i.e., "on-line". A heuristic criterion 
uses the on-line finite averages to terminate the trajec- 
tory integration once a desired tolerance on convergence is 
achieved. As a result, different initial conditions will result 
in different lengths of trajectories, depending on the regu- 
larity of trajectory that evolves from an initial condition. 

The ODE model of a system is numerically integrated 
producing state vectors ipt n (x) at adaptivcly-determincd 
time instances t n . Summing a zcroth-ordcr interpolation of 
evolution of fk approximates the continuous-time average 
along the trajectory: 
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To determine the stopping time of the integration, we 
pause the integrator in intervals of fixed length T e and at 
every time instance to + jT e compare the vector of time 
averages with the vector obtained at the previous pause 
to + U - l)T e : 



ADIFF J (x)= Fj(x) - F^x) 



where Fj (x) indicates the vector containing all the chosen 
obscrvablcs averaged along trajectory starting at x, com- 
puted in the [t ,t + jT e ] interval. When ADIFF^x) < 
ATOL for some pre-specified small tolerance ATOL, the 
integration is stopped, and Fj (x) taken to be the value 
of time-averaged observables at point x. The underlying 
assumption is that the finite-time averages will show little 
variation as the infinite limit is approached. As an alter- 
native to || - 1| the H~ s distance (J3|) could have been used 
in the stopping criterion, however, we have not explored 
that option here. 

Theory gives some guidance on times required for con- 
vergence of finite-time averages: over regular trajectories, 
averages converge with T -1 , while over mixing trajecto- 
ries, the slope is gentler, T" 1 / 2 . A classical result asserts 
that for other types of trajectories the slope T~ a can be 
arbitrarily small 32, §3.2B] . A computational study of con- 
vergence with various slopes is given in (25| . From a prac- 
tical standpoint, the outlook is not so bleak: slopes of 
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small a are to be typically expected to occur along tra- 
jectories that pass close to homoclinic and heteroclinic or- 
bits, and get entrained around periodic orbits that exist 
in their vicinity (intermittency, or stickiness) [3l| . Regions 
of intermittency are typically small, e.g., for a pendulum 
perturbed by a small periodic forcing of magnitude e, the 
area of intermittent regions scales as elog(l/e) |37j |. 

In Figure IA.9I we present several metrics which illus- 
trate how well the time averages converge at the time the 
simulation is stopped by the ADIFF(x) < ATOL criterion. 
Time average vectors Foe (x) and -Fatol (x) both represent 
finite-time averages, computed for trajectories originating 
at initial conditions x. The difference between them is in 
the moment of termination of the averaging process. The 
vector Fatol (x) was obtained by stopping the averaging 
at time T(x), dependent on the initial condition, when the 
described stopping criterion reached ATOL = 1CT 5 . The 
vector F aD (x) mimics the "true" infinity, where averaging 
was terminated at the initial-condition-indcpendcnt time 
= 2.5 x 10 4 . We compared these finite-time vectors in 
the L°° norm l-jl^ (Figure [A.9bp and in the H~ s norm 
|| - 1| 2 _ s © (Figure lA.9c[) . as the latter is the relevant error 
for the remainder of the computation. 

Generated images demonstrate that the described ter- 
mination criterion results in stopping times that vary by an 
order of magnitude (Figure l"A.9a[) . As a result, the errors 
in convergence are consistently low, for regular, chaotic 
and intermittent trajectories. If we compare the magni- 
tude of errors in H~ s - and L°°-norms, it is clear that the 
use of H~ s norm, which discounts higher wavenumbcrs, 
has a desired side-effect of increasing the robustness of the 
averaging to finite stopping times. 

Even though the proposed criterion is only one of many 
criteria which could be used for stopping, we believe that, 
based on the results presented, it provides a good heuris- 
tic that consistently produced low errors in time averages, 
while keeping the simulation times reasonable. 

Appendix B. A finite approximation to empirical 
metric 

A finite approximation to the norm on the Sobolev 
space H~ s ([3]) is computed by truncating the infinite sum 
over wavevectors k £ Z D to a box k £ [— K, K] D . In 
this appendix we show that such approximation converges 
with rate \f~K with the rate constant depending on the 
dimension of the system. 

Proposition 1 (Error in fmite-wavenumber approxima- 
tion to H~ s norm). Let the error in approximation Ek be 



given by a weighted sum of coefficients, Ck = 
summed outside the [— K, K] D box: 



fk{x) - f k (y) 
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where s = (D + l)/2, and < Ck < C(D) are differences 
in averages of Fourier observables scaled to unit Li norm. 
Then 
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where £(D) is the error decay rate constant. 

Proof. Norms on IP are ordered in scale, || A;|| 00 < ||fc|| 2 on 
7L D , enabling us to bound Ek by summing over oo-spheres 
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where a K = ff{k £ IP : ||fc|| < k} count the number of 
wavevectors in each oo-ball. To bound a K , we expand the 
polynomials in k: 
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as there is less than D/2 nonzero coefficients, all bounded 
by the largest coefficient 2 ( ^jy 2 j ) ■ 

Choosing order s = (D + l)/2, the error is bounded 
above by tails of sums of inverse squares 
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To bound the decay constant Q(J)), we use an upper 
bound for the binomial coefficient, computed in (35| : 
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for m, n,p £ N and m > p. Setting n = [D/2\ + 1, m = 2, 
p = 1, provides a conservative estimate: 

/ D \ [8 2 D 
\[D/2\) < V 7^/DTT■ 

As observables are /^-normalized Fourier functions (|4|), 
C(D) = (2tt)~ d / 2 , allowing us to compute a bound on 
£\D): 
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Figure A. 9: Convergence of averages for N = 500 trajectories of ABC flow starting at plane z = 0. Simulations were 
terminated at tolerance ATOL = 1 x 10~ 5 . Figures show: IA.9al the termination time Tlx). lA~9bl IA.9cl the difference 
between time averages at time T(x) and at "true" infinity = 2.5 x 10 4 , in infinity and H~ s norms, respectively. 
Obscrvables were Fourier harmonics d3| with wavevectors k e [—3, 3] 3 , and simulation extension time was T e = 500. 
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Appendix C. Diffusion Maps 

The Diffusion Maps is an algorithm that retrieves an 
orthogonal coordinate system, the diffusion coordinates, 
for a set of points in a metric space. If the set of points 
is taken from a differentiable submanifold in the original 
space, the retrieved coordinates are discretized eigenfunc- 
tions of a Laplacc-Beltrami operator on the submanifold. 
Diffusion coordinates yield a parsimonious description of 
the data set, i.e., even a truncated diffusion coordinate set 
provides a good low-dimensional model of the original data 
set. 

The Diffusion Maps is a part of a broader class of lapla- 
cian manifold-learning algorithms, which rely on spectral 
computations to obtain information about geometry of 
manifolds. The algorithm was first described by Lafon 
and Coifman in [(1,13] with more detailed discussions on 
numerical implementation given in 23J, |24| . To convey the 



intuition about the algorithm, we provide a short account 
of the general algorithm and implementation details we 
found relevant. At the end, we discuss how we applied it 
to the ergodic quotient. 

The ambient space is the space C N with a distance 
function T>. We aim to obtain a representation of a subset 
y C <C N , with a distribution dP(y) = p(y)dy supported on 
y. Averaging operators Ah are defined through the kernel 
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and they represent diffusion of points on y, while the dif- 
fusion time-scale h specifies the width of the gaussian ker- 
nelH The h — > limit of the flow A* h has a simple struc- 
ture when y has a well-defined differential structure: it is 
generated by a Fokkcr-Planck-type operator G: 
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where diffusion evolves over time period t. The Diffusion 
Maps algorithm computes values of eigenfunctions (y) 
of e~ tG at a finite number of points j/j that are sampled 
from y according to the (sampling) density p(y). 

Eigenfunctions of the Laplacc-Beltrami operator con- 
vey the geometry of the set y. If understanding geome- 
try of y is the goal, as in the case of ergodic quotient, it 
is necessary to remove the bias in G resulting from non- 
uniformity in density p. The bias can be removed by re- 
placing the heat kernel kh by its rescalcd version: 



kh 



p(y) 



kh(y,z) 

p(y)p(z) 



kh(y,z)dz, 



where p acts as an estimate for the sampling density. In 
this case, the flow in the h — > limit represents the pure 
Laplacc-Beltrami diffusion on y, i.e., G = A. 



5 Sometimes 1/h is referred to as the bandwidth. 
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To compute the diffusion coordinates, denned as the 
embedding of points yi into eigenvalue-scaled eigcnfunc- 
tions of the diffusion operator y — > {Xix'^iv), X2X^(y), ■ ■ ■), 
Diffusion Maps treats the finite set of sampling points 
j/j G C N as vertices of a fully-connected graph in y, while 
weights of edges are given by the ambient distance T> on 
C D , between pairs of points j/j. The diffusion process on 
the graph can be represented by a Markov chain, whose 
eigenvectors x^ are discretizations of eigcnfunctions of 
the diffusion on y. In other words, we construct a discrete 
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From the standpoint of the implementation, Diffusion 
Maps takes a matrix of pairwise distances between points 
yi as the input, and computes diffusion coordinates A% 
following these steps: 

V matrix (<£) = V 2 {y u y 3 ), (C.3a) 



Heat kernel Ay = exp 
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Unbiased heat kernel A,- 
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Solve the eigenproblem S\ = Ax, 



(C.3b) 

(C.3c) 

(C.3d) 

(C.3e) 
(C.3f) 



with i, j € 1,2, . . . , JV indexing elements of matrices and 
vectors, and T>- s _k representing summation ([3]) truncated 
to k G [-K, K] D . 

The time-scale of the heat kernel h in (|C.3b|) is a pa- 
rameter that determines the strength of diffusion along 
graph that samples y. It influences the information about 
topology of y that we infer from discrete data: too small 
of a choice of h will introduce artificial discontinuities, 
while too large of a choice will artificially connect dis- 
connected components of y. Paper [24], §5.3] describes 
several heuristics for determining h from the data set; we 
chose the Neighborhood Size Stability approach. It com- 
putes the minimal time-scale h for which every graph ver- 
tex has N m i n neighbors within the characteristic diffusion 
distance \/2h, as measured by the H~ s distance. Num- 
ber N m i n is chosen by the user, depending on number of 
discrete samples analyzed. 

Since S is not symmetric, complex values appear as 
partial results in numerical solutions of its eigenproblem, 
which is a source of further numerical errors. To avoid 
such issues, it is common to first symmetrize the matrix, 
using a spectrum-preserving symmetrization: 



A* 



Solving the eigenproblem for S, we obtain eigenvectors 
yX k ) from which we recover the eigenvectors that sam- 
ple the diffusion eigenfunctions (|C.3fp by rescaling x^ by 
the zeroth eigenvector x\ = Xi/xi - Vectors x^ are 
point- wise evaluations of diffusion eigenfunctions x (y) 

at yi, i.e., Xi = X (k) {Vi)- 

In application of Diffusion Maps to the ergodic quo- 
tient, initial conditions x% G Ai and associated trajectories 
get mapped to the points yi in the space C N of averaged 
observables, with dimension N given by the number of 
chosen observables. The ambient distance T> is the H~ s 
distance, while the set y is the ergodic quotient £. The 
sampling density p(y) is determined both by the distri- 
bution of initial conditions Xi of trajectories on the state 
space, which can be controlled, and the distribution of the 
ergodic sets in the state space, which is a priori unknown 
and cannot be controlled. As a consequence, p(y) is rarely 
uniform and the rescaling of the diffusion kernel is a nec- 
essary step in the algorithm to obtain the intrinsic coor- 
dinate set with the correct interpretation of the euclidcan 
distance. 
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